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Abstract 


We present a new method to diagnose the middle atmosphere climate sensitivity by 
extending the Climate Feedback-Response Analysis Method (CFRAM) for the coupled 
atmosphere-surface system to the middle atmosphere. The Middle atmosphere CFRAM 
(MCFRAM) is built on the atmospheric energy equation per unit mass with radiative heating and 
cooling rates as its major thermal energy sources. MCFRAM preserves the CFRAM unique 
feature of an additive property for which the sum of all partial temperature changes due to 
variations in external forcing and feedback processes equals the observed temperature change. In 
addition, MCFRAM establishes a physical relationship of radiative damping between the energy 
perturbations associated with various feedback processes and temperature perturbations 
associated with thermal responses. MCFRAM is applied to both measurements and model output 
fields to diagnose the middle atmosphere climate sensitivity. It is found that the largest 
component of the middle atmosphere temperature response to the 11 -year solar cycle (solar 
maximum vs. solar minimum) is directly from the partial temperature change due to the variation 
of the input solar flux. Increasing CO 2 always cools the middle atmosphere with time whereas 
partial temperature change due to O3 variation could be either positive or negative. The partial 
temperature changes due to different feedbacks show distinctly different spatial patterns. The 
thermally driven globally averaged partial temperature change due to all radiative processes is 
approximately equal to the observed temperature change, ranging from -0.5 K near 25 km to 
-1 .0 K near 70 km from the near solar maximum to the solar minimum. 
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1. Introduction 


The warming of Earth’s surface and lower atmosphere is associated with enhanced 
middle atmosphere cooling and a strengthening of the Brewer-Dobson circulation through 
radiative, dynamical, and photochemical coupling. Because both the air density and the optical 
depths of major radiatively active species decrease with altitude, the physical state of the middle 
atmosphere as represented by various parameters such as temperature and winds is quite 
sensitive to climate forcing and is thus a good indicator of surface global wanning. Hence, a 
more accurate quantification of the middle atmosphere responses to solar variability and 
anthropogenic changes in trace species is necessary to improve predictions of climate change. 

The Climate Feedback-Response Analysis Method (CFRAM) has been developed for 
separating and estimating various climate feedbacks in the coupled troposphere-ocean system 
(Fu and Cai 2009; Cai and Fu 2009; hereafter FC09 and CF09). CFRAM is formulated based on 
the atmosphere-surface energy equation, and it explicitly decomposes the directly measurable 
total temperature change into partial temperature changes due to individual external forcing and 
feedback processes (FC09, CF09). The unique feature of CFRAM is that this decomposition into 
partial temperature changes is locally additive, so that the total temperature change is the sum of 
all the partial temperature changes at every spatial point. From the modeling perspective, the so- 
called external forcing and its variation of a system are akin to independent variables or 
parameters that would be prescribed as input values in a model. On the other hand, the feedback 
or internal processes of a system are similar to dependent variables or parameters that often 
constitute a set of model output values. 

In this paper, CFRAM is extended to the middle atmosphere based on three physical 
features of this region: (i) radiative energy exchange plays a major role in the energy budget; (ii) 
the air density varies with altitude by several orders of magnitude and the energy deposition per 
unit mass is often scaled by a factor that slowly varies with altitude or log-pressure; and (iii) the 
energy flux associated with the level of the Earth’s surface and the layered middle atmosphere 
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are not directly coupled. As a result, the Middle atmosphere Climate Feedback-Response 
Analysis Method (MCFRAM) is formulated by the energy equation in a fonn of heating and 
cooling rates per unit mass in a commonly used unit of K day 1 . Its mathematical form is similar 
to a well-documented radiative transfer technique for analyzing radiative damping or relaxation 
of the atmospheric temperature disturbances ( e.g ., Goody and Yung 1989, Zhu and Strobel 
1991). The newly developed MCFRAM is here applied to the middle atmosphere to derive 
various partial temperature changes based on both satellite measurements and output of a three- 
dimensional (3D) chemistry-climate model (CCM). 

In Section 2, we briefly review and extend CFRAM to the middle atmosphere. Then, we 
perfonn the fundamental eigenmode analysis to the generalized damping matrix derived from the 
MCFRAM. The middle atmosphere temperature and ozone fields needed in the analysis are 
derived from the Sounding of the Atmosphere using Broadband Emission Radiometer (SABER) 
onboard the Thennosphere, Ionosphere Mesosphere, Energetics and Dynamics (TIMED) 
satellite. Section 3 shows the MCFRAM results derived from the SABER measurements whereas 
section 4 performs a set of similar MCFRAM analyses on the output fields of the Goddard Earth 
Observing System chemistry-climate model (GEOSCCM; Pawson et ai, 2008, and references 
therein). Section 5 summarizes the paper. 

2. Review and extension of the Coupled Feedback Response Analysis Method 


2.1 Formulation of the middle atmosphere CFRAM 

CFRAM was originally formulated in a form of a vertical energy flux difference for a 
single-column energy equation in a form of the time mean energy balance equation (LC09; 
CL09): 




( 1 ) 


where R and S ' are the infrared and solar flux differences corresponding to total radiative 
cooling and heating of a layered atmosphere, respectively. Q is the non-radiative energy flux 
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convergence in the atmospheric layers. T is temperature profile, ( r , s, ...) are the mixing ratios of 
radiatively active species such as CO 2 , O 3 , H 2 O and clouds, and (a, /?, ...) are the parameters 
such as the solar irradiance, surface albedo and solar declination angle that will affect the 
atmospheric energy. The terms in the energy Eq. (1) for CFRAM have the units of energy flux 
W m 2 , which corresponds to the heating or the cooling rate per unit volume for a given layer of 
atmosphere. There are several advantages of adopting the flux form with units W m 2 in the 
classic CFRAM: (i) the energy flux of the atmosphere can be naturally coupled with the surface 
energy flux; (ii) the top of the atmosphere (TOA) version of CFRAM can be directly compared 
to a TOA-based climate feedback analysis such as the partial radiative perturbation (PRP) 
method; (iii) the layer thickness of the tropospheric CCMs is usually slowly varying in mass so 
the heating or cooling rate perturbations per unit space of different layers also slowly vary with 
altitude. 


The air density decreases with altitude exponentially in the middle atmosphere, ranging 
from the tropopause (~10 km) to the turbopause (-110 km), spanning several orders of 
magnitude in density variation. The energy deposition or the atmospheric heating rate in 
measurements and models is often scaled in a unit mass with a setting in vertical grid that slowly 
varies with altitude or log-pressure. As a result, we begin by developing our MCFRAM from an 
energy equation per unit mass, i.e., by dividing Eq. (1) by c ptsz with c p , p and Az being the 
specific heat at constant pressure, air density and layer thickness, respectively, 


R(T,r,5,...,a,y0,...) = S(T,/y5,...,«,y0,...) + Q(T,r,5,...,«,y0,...) + Q mo/ (T), 


( 2 ) 


where R and S are the radiative cooling and heating rates, respectively. Q is the non-radiative 
heating rate excluding the molecular thermal conductivity Q mo/ (T) that is only a function of 
temperature profile T (Banks and Kockarts 1973). The units of all terms in Eq. (2) are K day -1 . 
We now consider two statistical equilibrium states 1 and 2 with two different sets of 
corresponding atmospheric parameters all satisfying the energy balance equation (2). In practice, 
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these two states can be two ensemble, time or spatially averaged mean states. The difference of 
the energy equations between these two states is 

A(R-Q fflo/ ) = AS + AQ. (3) 


We now introduce the linear approximation to the responses of R and Q mol to the temperature 
variation and separate this term from the variations due to other parameters: 


A(R-Q,„ o/ )* a(R 3 o/) AT + [R(T,r 2 ,s 2 , /? 2 , s lf /?„ -)] 

<3T 


(4) 


where T is the “mean temperature profile” between profiles Tj and T 2 . Substituting Eq. (4) into 
Eq. (3), we obtain 


AT = A -1 {AS - A'R + AQ} , 


(5) 


where A = 3(R- Q mol ) / 0T « A(R - Q mol ) / AT is the generalized damping matrix in units of 
day -1 and A'R defined by the last two terms in Eq. (4) is the change in total cooling rate due to 
all parameters except the temperature profile. 

In CFRAM, where the surface and the atmosphere are strongly coupled radiatively and 
dynamically, the discretization of the energy equation (1) and the derivation of the “Planck 
feedback matrix” 0R 5T based on the temperature profile need to include the temperatures of 
both the surface level and layers of the atmosphere (LC09). The surface temperature and 
atmospheric temperature are treated as equally important in the setting of the problem. As a 
result, Eq. (1) needs to be in the form of energy flux divergence in units of energy flux W m 2 . 
The temperature profile in the middle atmosphere is not directly coupled to the Earth surface. It 
can therefore be discretized solely based on a layered atmosphere in an energy equation of unit 
mass (Eq. 2). The effect of the energy flux emergent from the lower boundary on the middle 
atmosphere is primarily the radiative flux that is independent of the temperature in the middle 
atmosphere. For example, the effect of the solar radiative flux can often be parameterized by an 
effective albedo of the surface and lower atmosphere ( co 0 ), which enhances the heating rate due 
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to absorption of the Chappuis bands (410-750 nm) by ozone in the stratosphere caused by 
surface reflection and multiple scattering of clouds, aerosols and air ( e.g ., Meier et al. 1982; 
Nicolet et al. 1982). The calculation of the generalized damping matrix A for a basic state of 
temperature and species distributions can be implemented by a radiation algorithm and molecular 
diffusive fonnulation. In this paper, the Johns Hopkins University Applied Physics Laboratory 
(JHU/APL) middle atmosphere radiation algorithm (Zhu 1994, 2004) is adopted for radiative 
cooling calculations, and a temperature-dependent thermal conductivity of A = 5.6x1 0^T ° 69 
[kg-m-s -3 -KT 1 ] (Banks and Kockarts 1973) is used for calculating the diffusive heat flux of 
AST / dz . Each column (vertical axis) of A represents the vertical profile of cooling rate and 
diffusive heating rate difference ( K day 1 ) due to a unit change in temperature at altitude z 
(horizontal axis). 

In the middle atmosphere, the effect of line overlap is negligible for the infrared radiative 
cooling rate calculations. As a result, the total infrared cooling rate can be evaluated as the sum 
of the cooling rates due to CCE, O3 and H2O (Zhu 1994). Therefore, the term A'R in Eq. (4) or 
(5) becomes, 

A'R = AR C02 + AR 03 + AR H2 ° . (6) 

Note that for the middle atmosphere, Eq. (6) is nearly exactly satisfied. In other words, the linear 
separation of the partial infrared radiative cooling rate due to individual gases in the middle 
atmosphere is satisfied automatically, which is not the case for the troposphere. Therefore, the 
only linear approximation introduced to the infrared radiative cooling rate is the separation of the 
temperature variation, as indicated in Eq. (4). In this sense, fewer approximations for the middle 
atmosphere feedback analysis have been used than those for the troposphere and surface 
temperatures, e.g., the CFRAM (LC09) and PRP method (Soden et al. 2008). 

For radiative heating by solar flux, we still need to invoke a linear approximation to 
decompose the energy perturbation into individual components by different factors, namely, 
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ac 3S 6S i^S 3S 

AS ~ _ — — A[0 3 ] + — — A[O 0 ] + — — — A/^jgy + — — A co. 


S[0 3 ] 


d[0 2 ] 


dR 


107 


dco n 


0 ’ 


(7a) 


or 

AS « AS 03 + AS 02 + AS f107 + AS'" 0 , (7b) 

where [O 3 ] and [O 2 ] are the ozone and oxygen densities, respectively. F m is the 10.7-cm solar 
radio flux (in units of 1 0 22 W m 2 Hz 1 ), which is a parameter representing the solar flux 
variations, and co 0 is the effective albedo of surface and the lower atmosphere. The effect of [O 2 ] 
variation on the energy perturbation (AS 02 ) is only important in the lower thermosphere. In the 
middle atmosphere, the vertical velocity associated with the meridional circulation or the 
residual circulation plays an important role in coupling the radiation with dynamics and 
photochemistry. Here, invoking a linear approximation, we may explicitly extract this special 
term of “dynamical response” from the total non-radiative energy source (Holton 2004): 

AQ = 0Aw* + (A0)w* + A"Q eddy , (8) 

where the diagonal matrix 0 is the static stability parameter and column vector Aw is the 
variation in resolved vertical velocity that yields the change in adiabatic cooling. The last term in 
Eq. (8) represents the contributions due to unresolved small-scale eddies and the energy transport 
by horizontal wind among neighboring vertical columns. Note that although the non-radiative 
energy source ( AQ ) can be evaluated from the dynamical modules during runtime of model 
integrations as reported in Lu and Cai (2010) and Song et al. (2014), it cannot be obtained 
directly from observations. It can also be evaluated as an energy residual term to balance the net 
radiative cooling rate and molecular mixing according to Eq. (3): AQ = A(R - Q mo/ - S) . Such an 
approach of using better-defined thermal forcing to diagnose mechanical forcing was also 
proposed in Zhu et al. (2001) to diagnose the dynamical fields in the middle atmosphere. As 
reported in Lu and Cai (2010) and Song et al. (2014), AQ inferred explicitly from dynamical 
fields is almost identical to that inferred as an energy residual term. Given AQ , we then use Eq. 
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(8) to obtain A"Q <Y/ ' /l from the difference between AQ and the other two terms which can be 
calculated from the available T and w profiles. 

Substituting Eqs. (6)-(8) into Eq. (5), we obtain 

AT = Z • {-AR C02 + (AS 03 - AR 03 ) + AS F107 - AR H2 ° 

+AS" 0 + 0Aw* + (A0)w* + A"Q eddy } , (9) 

where Z = A -1 is the generalized relaxation matrix. Note that the change of [O3] in the middle 
atmosphere contributes to both radiative heating and cooling rate variations. This is similar to 
H2O and clouds in the troposphere that can both radiatively heat and cool the atmosphere. 

As indicated in LC09, Eq. (9) represents the property of additive thermal responses, i.e., 
the sum of the partial temperature changes ( AT (n) ) of the MCFRAM due to individual variations 
of various external forcing such as A[C0 2 ] or AF 107 and feedback processes such as A[0 3 ] is 
the total temperature change ( AT' oto/ ): 

AT' oto/ = ^ AT ("), (10) 

n 

where 

AT ( " ) =Z-AF'" ) . (11) 

Here, we define the total temperature change AT' 0 '" 7 to be an observed quantity representing the 
actual difference in the measured temperatures between the two equilibrium states. The energy 
perturbations AF (,!) in Eq. (1 1) denote various terms in the brackets on the right-hand side of Eq. 

(9) . Physically, AT <H) correspond to the partial temperature changes associated with linear 
atmospheric thermal responses to the energy perturbations AF (n) caused by individual parameter 
variations. Those parameter variations can either be derived from observations or from model 
output. The physical meanings of these partial temperature changes (AT <n) ) in MCFRAM are 
shown in Table 1. The sum of the first six components forms the partial temperature change due 
to radiative processes AT""' = AT (I 6) . The non-radiative partial temperature change AT""” rad 
includes changes due to both the grid resolved and unresolved atmospheric motions. It should be 
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noted that AT'™' has been derived from the changes in net radiative heating rate excluding the 
cooling rate change due to the temperature variation itself, i.e., the terms (AS-A'R ) in Eq. (5). 
One important reason that the temperature variation is singled out in Eq. (4) is that the 
generalized damping matrix A introduced in Eq. (5) is well-behaved and always invertible. 

The additive relation (10) for the temperature changes is an alternative expression of the 
energy Eq. (3) that too is additive. A linear transformation that singles out the total temperature 
change from energy difference on the left-hand side of Eq. (3) also leads to the partial 
temperature differences as shown in Eq. (11) and allows us to derive this alternative relationship. 
The principal advantage of the additive relation (10) for temperature over the additive relation 
(3) for energy is that AT'"'" 7 on the left-hand side of Eq. (10) is a directly observed and 
commonly used quantity, which can serve as a natural and standard scale for comparison. 

In addition to temperature T and its changes AT <n) , we may also choose a different 
variable to express the energy budget and its changes. For example, given a vertical profile of 
“thermal response” as expressed by the temperature changes AT""" / and AT ( ” ) , we may calculate 
the corresponding local “dynamical response” of changes in meridional circulation by the 
following linear transformation 


Aw' ota/ =(0“ 1 A)-AT 


total 


and 


Aw (,!) = (0 _1 A)-AT' 


O) 


( 12 ) 


(13) 


where we have already assumed a stable stratification of the atmosphere so that the matrix 0 
never becomes singular or ill-conditioned (Holton 2004). This condition generally holds well in 
the middle atmosphere. Substituting Eqs. (12)-(13) into Eq. (10) yields a different form of energy 
budget equation: 


Aw tota> = Yj Aw (,l) . 


(14) 
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Again, note that in MCFRAM the total temperature change AT'"' a/ is a physical quantity that is 
directly measurable. On the other hand, the total resolved local vertical velocity defined by Eq. 
(12) is only an equivalent quantity corresponding to the observed AT"" o/ . The physically 
measurable vertical velocity is Aw" and is related to the contribution of an equivalent partial 
temperature change via Aw =(0 A)-AT :: (Table 1). 

It is noted that the global “thermal response” and “dynamical response” are closely 
related in a more general perspective under the statistical equilibrium condition (Fels 1987; Zhu 
et al. 2001). The meteorological underpinning of such a relation in a meridional plane is the 
thermal wind balance. Specifically, even though the meridional circulation is driven by the 
meridional gradient of the diabatic heating, the vertical gradient of the diabatic heating is exactly 
balanced by the meridional gradient of the mechanical forcing (Fels 1987; Zhu et al. 2001). As a 
result, the strengthening of the meridional circulation such as the Brewer-Dobson circulation in 
the lower stratosphere can be interpreted either as a response to changes in thermal forcing or as 
a response to changes in wave drag. Therefore, the spatial structures of AT 1 ' 0 derived from 
MCFRAM based on the energy equation also provide us with a global insight into both the 
thermal and dynamical responses in the middle atmosphere. Note that MCFRAM as outlined by 
Eqs. (10)-(11) together with Table 1 generally applies to independent columns of the middle 
atmosphere. The spatial structure of the derived AT (n) is only a result of the diagnostic analysis 
but is not explicitly included in the analysis procedure. 

2.2 Eigenmodes of the generalized damping matrix and illustration of MCFRAM 

In Eq. (5) or (9), there is a common matrix factor that linearly multiplies all the radiative 
and non-radiative energy perturbation terms. As a result, both the magnitude and vertical 
structure of the climate feedbacks are significantly influenced by the generalized damping matrix 
A defined in Eq. (5) or the generalized relaxation matrix Z defined in Eq. (9). Table 1 explicitly 
shows that the partial temperature changes are proportional to the energy perturbation vectors 
AF < " ) for different processes and are modified by the same generalized relaxation matrix Z. 
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Specifically, for a given vertical profile of the energy perturbation the spatial structure of partial 
temperature change is completely determined and can be understood by the eigenvectors ( % l ) 
and eigenvalues ( A,. ) of A or Z = A 1 : 


A $ i =A£ l or Z^=X~%, 1=1,2, 


(15) 


where N is the total number of vertical layers. Equation (15) indicates that the eigenvalues of Z 
are the inverse of the eigenvalues of A corresponding to the same eigenvectors. Here, Z ( . and 2 1 
can be called generalized damping rate and relaxation time corresponding to the perturbation 
eigenvector , respectively. In the absence of molecular viscosity ( Q mol = 0 ) the generalized 
damping matrix is given by A = SR / ST . Its eigenvalue 2,. happens to be the radiative damping 
rate of a temperature perturbation ( e.g ., Goody and Yung 1989, Zhu and Strobel 1991). The 
effect of the vertical structure of the temperature perturbation characterized by its eigenvector 
on the radiative damping rate has been well documented (Zhu and Strobel 1991, Zhu 1993). The 
occurrence of the radiative damping rate in MCFRAM is a natural consequence that the basic 
MCFRAM equation (9) or (10) is an energy perturbation equation. When the energy perturbation 
is specifically referring to the cooling rate change in association with a temperature perturbation 
that has been singled out among all the other changes, it is the radiative damping rate that 
establishes the connection between the two perturbations. In general, the magnitude of 2 ( . under 
non-vanishing Q mo/ conditions is proportional to the magnitudes of the radiative cooling rate and 
molecular viscosity. It increases with the increasing characteristic vertical wavenumber of the 
energy perturbation, i.e., the wavenumber of cooling rate variation or the temperature variation. 

The infrared radiative heat exchange by CO2 and O3 makes a major contribution whereas 
cool-to-space cooling by H 2 O makes a minor contribution to the radiative cooling rate in the 
middle atmosphere (Zhu 1994). Here, we use the T and [O 3 ] observed from the SABER onboard 
the TIMED satellite to derive A or Z and to perform the eigenmode analysis to illustrate the 
general characteristics of the eigenvector of A or Z in the middle atmosphere. The needed global 
mean H 2 O profile for the radiation algorithm is derived from the 3D Goddard Earth Observing 
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System chemistry-climate model (GEOSCCM; Pawson et al., 2008). In Fig. 1, we show the 
TIMED/SABER measured global mean T and [O3] averaged over a 54°S-54°N latitudinal range 
and a 12-year period of 2002-2013. The SABER measurements ranging from 20 km to 110 km 
in the middle atmosphere are merged with the US Standard Atmosphere in the troposphere. The 
vertical resolution of all the input profiles from surface to 1 10 km is about 0.7 km. The radiative 
heating and cooling rate calculations based on the JHU/APL radiation algorithm are performed 
in the entire vertical domain of 157 layers whereas the MCFRAM analysis is applied to the top 
129 layers (7V=129) that corresponds to a middle atmosphere ranging from 10 km to 110 km. 
The matrix A or Z has dimensions of 129 x 129 with 129 eigemnodes. Any given vertical profiles 
of the energy perturbations (AF (,!) ) can be decomposed by a complete set of the eigenvectors, 
with each component decaying, i.e., relaxing to 0, at a rate proportional to the inverse of their 
corresponding eigenvalues. Figures 2 shows a set of 17 selected vertical eigenmodes of the 
generalized damping matrix A calculated from T and [O3] shown in Fig. 1 based on the 
JHU/APL radiation algorithm (Zhu 1994, 2004). The CO 2 volume mixing ratio in the calculation 
is set at a 2005 level of 380 ppmv. The eigenmodes describe a quantitative relationship between 
the energy perturbations and the corresponding temperature perturbations. The eigenvalues ( X t ) 
of the selected eigenvectors (2j ; .) range from a maximum value of Z max =18.98 day -1 (blue line 
marked with circles in Fig. 2a) to a minimum value of A mm = 0.023 day -1 (blue dashed line in 
Fig. 2d). The vertical eigenmode of the largest damping rate corresponding to the smallest 
relaxation time (X maK ~ l =0.053 day) is a wave packet located at 107 km with a very small 
vertical scale of ~4 km. That particular wavy energy perturbation will be effectively smoothed 
out in a very short period and produces a very small temperature perturbation. On the other hand, 
the eigenmode with the smallest damping rate has a vertical structure of a near uniform heating 
or cooling near the tropopause. This mode has the largest relaxation time (X mm ] = 43.4 day) that 
will yield the largest response in partial temperature change for a given unit of heating or cooling 
rate perturbations. 
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There are two distinct features shown in Fig. 2. First, there exists a strong scale- 
dependence of the eigenvectors for the generalized damping matrix A. Eigenvectors 
corresponding to large-scale vertical perturbations have small eigenvalues. Second, the 
magnitude of the eigenvalue decreases as the location of the characteristic perturbation shifts 
from the upper middle atmosphere to the lower middle atmosphere. As a result, we note that 
when the value of eigenvalue decreases as we move consecutively from panel (a) to panel (d) the 
vertical scale of eigenvector increases and the location of its main perturbation shifts to the lower 
altitude. This is consistent with the general nature of the radiative damping of temperature 
perturbation in the middle atmosphere ( e.g ., Goody and Yung 1989; Zhu 1993). In addition, the 
effect of the molecular diffusion included in A has the same general characteristics of small-scale 
perturbations at a higher altitude being more effectively damped or filtered. To show the general 
nature of scale-dependence and its departure from a precise one for the eigenmodes in the entire 
middle atmosphere we perform a Fourier transform to all 129 eigenvectors and calculate their 
power spectral densities (PSDs) (Zhu and Strobel 1991; Zhu 1991). Figure 3 shows a scatter plot 
between the generalized damping rate A ; and the wavenumber of the maximum peak in the PSD 
for all 129 eigenvectors. Also shown in the figure are the analytic expression for the 
parameterized radiative damping rate proposed in Zhu (1993) and a square fit (A,. ~ nr) to the 
diffusive damping. Because of the vertical inhomogeneity of the atmosphere, the relationship is 
not single-valued. For example, a wave packet with a large vertical scale located in the 
mesopause could have the same damping rate as one in the stratosphere with a small vertical 
scale. A better parameterization for radiative damping in practice is to introduce a scale- 
dependent radiative damping rate that also varies with altitude (Fels 1982; Zhu 1993). 

The effect of the vertical structure of the energy perturbations on the partial temperature 
changes through A or Z can be seen from Fig. 4 where the partial temperature changes of AT C02 , 
AT 0 ' and AT FI07 are calculated from Table 1 based on three energy perturbations caused by 
changing three atmospheric parameters (i) CCF volume mixing ratio is doubled from 380 ppmv 
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to 760 ppmv, (ii) O 3 volume mixing ratio is uniformly reduced by 50%, and (iii) solar index F 10.7 
is increased from 60 to 260 (in units of 10 22 W m 2 Hz ). The vertical structure of temperature 
difference (Fig. 4b) is smoother than and significantly different from that of the heating rate 
variations (Fig. 4a). This is mainly due to the scale-dependence of the generalized damping rate 
( ) where smaller scale energy perturbations are more effectively damped, i.e., partial 
temperature changes are smoother than energy perturbations. Furthennore, the lower middle 
atmosphere is more sensitive in partial temperature changes to a smaller energy perturbation due 
to greater opacity than the upper middle atmosphere. 

3. Application of MCFRAM to TIMED/SABER measurements 

Application of MCFRAM is straightforward once the input fields of various parameter 
variations as indicated in Eqs. (5) and (9) together with Table 1 are available. While climate 
models (such as the GEOSCCM) can provide all the needed and uniformly distributed global 
input fields, satellite measurements often provide only part of the needed fields to derive the 
balanced additive relation (10). In this section, we show the MCFRAM analyzed results by using 
SABER measured T and O 3 fields (Russell et al. 1999). We use the VI. 07 SABER data available 
to the public from the TIMED mission data center (http://www.timed.jhuapl.edu) which yields 
significantly improved temperature retrievals at high latitude summer (Kutepov et al. 2006). 

Figures 5a and 5b show the zonal mean T and O 3 fields in the middle atmosphere derived 
from SABER measurements in the low and mid-latitudes averaged over a 12-year period of 
2002-2013. Shown in Figs. 5c and 5d are the T and O 3 difference between two time-mean states 
covering the periods of 2002-2003 and 2008-2009, respectively. Though the overall temperature 
in the middle atmosphere exhibits a noticeable decrease from the 2002-2003 period of near solar 
maximum to the 2008-2009 period of solar minimum over most regions, there are some regions 
showing positive temperature anomalies in response to the solar energy input decrease. We note 
that the observed temperature difference represents the total effects contributed by various 
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processes including the solar flux changes due to solar cycle and man-made variations in CO2 
concentration and other chemical species. 

We now apply MCFRAM to the SABER observed T and O3 difference between two 
periods of 2002-2003 and 2008-2009. The corresponding mean CO2 mixing ratios and solar flux 
indices used in MCFRAM analysis for these two periods are [r C 02~374.7 ppmv, Fi 07 ~167.1] and 
[i'C 02 ~ 386.3 ppmv, F 107—68. 1], respectively. There are 12 yaw cycles in each two-year period 
with each yaw cycle covering about 60 days. The corresponding local time and latitudinal 
coverage in two yaw cycles separated by six years are nearly identical. The temperature 
difference shown in Fig. 5c represents the total temperature difference AT ,oto/ as defined in Eq. 
(10). The MCFRAM analysis is perfonned separately to the corresponding yaw cycles with the 
seasonal parameters such as the solar declination angle and F10.7 varying with different yaw 
cycles. The partial temperature changes as defined in Eq. (11) or Table 1 will be the 12-yaw 
cycle mean of all partial temperature changes between the two yaw cycles in the corresponding 
time periods separated by 6 years. Given the observed T, O3 and F10.7 variations and using the 
JHU/APL middle atmosphere radiation algorithm, the first three components of the partial 
temperature changes shown in Table 1, i.e., AT to2 , AT 03 and AT FI07 , can be explicitly 
evaluated. Since H2O and other radiatively active species only make minor contributions to the 
radiative heating and cooling rate in the middle atmosphere, we expect the sum of the above 
three terms is approximately the partial temperature change due to radiative transfer AT™ f/ as 
described in Table 1. As mentioned before, we use the energy residual of Eq. (3) to estimate AQ 
to calculate AT"°" _rai . In Fig. 6, we show the latitude-altitude distributions of AT C02 , AT 03 , 
AT FI07 and AT"°" _rad . Also shown in the figure are Aw C02 defined by Eq. (13) and the error in 
At non-rad 1 j ncar j za tj on , / e the difference between AT" 0 " rad based on the energy residual 

(Table 1) and the one based on a temperature residual AT' ote/ -(AT C02 + AT 03 + AT FI07 ) . 

We note that the middle atmosphere cooling rate by the CO2 1 5 -pm band is mainly 
contributed from its cool-to-space component with its escape probability slowly varying with 
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altitude in the middle atmosphere (Zhu et al. 1992). A uniform change in CO 2 mixing ratio also 
leads to a near uniform change in escape probability in the middle atmosphere. Hence, the 
maximum response to a uniform increase in CO 2 mixing ratio in the middle atmosphere occurs 
at the equatorial stratopause (Fig. 6 a), where the peak temperature as shown in Fig. 5a produces 
the largest cooling rate variation. On the other hand, the response AT 03 due to change in O 3 
concentration represents a combined effect of both the solar flux heating and O 3 9.6 pm band 
infrared cooling. Since there are both positive and negative ozone variations between 2002-2003 
and 2008-2009 periods (Fig. 5d), the induced partial temperature change AT 03 also shows a non- 
uniform spatial pattern (Fig. 6 b). The peak variation in temperature in the upper mesosphere is 
mainly due to the change in localized absorption of solar ultraviolet (UV) flux heating whereas 
the peak variations in the stratosphere are mainly due to the enhanced O 3 9.6 pm band cool-to- 
space cooling rate variations in a more transparent atmosphere. Here, we note that the middle 
atmosphere climate responses to the cooling rate changes induced by CO 2 and O 3 variations are 
different. AT C02 due to CO 2 variation (Fig. 6 a) mostly follows the total temperature field due to 
a strong dependence of outgoing infrared radiation on the Planck blackbody emission whereas 
AT 0 ' due to O 3 variation (Fig. 6 b) mostly follows O 3 concentration due to a stronger 
dependence of radiative emission on more rapidly varying escape probability (Zhu et al. 1991). 
AT FI07 shown in Fig. 6 c exhibits a pattern of overall monotonic increase in magnitude with 
altitude mainly due to the fact that solar UV fluxes of greater variations at shorter wavelengths 
are generally absorbed at higher altitudes. 

We note that the overall spatial pattern and magnitude of AT""" rad shown in Fig. 6 d is 
similar to AT """ 7 shown in Fig. 5c, indicating the importance of dynamical drive of the zonal 
mean middle atmospheric thermal structure. One striking feature in Fig. 6 d is that AT""" rad is 
significantly greater and has richer spatial structure than any individual partial temperature 
change due to radiation processes. In other words, most part of temperature changes in the 
middle atmosphere are associated with dynamic processes and the corresponding changes in 
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thermal radiation in turn balance the non-radiative energy source. One plausible explanation is 
that the middle atmosphere thermal radiative forcing as a whole is largely modulated by the 
dynamical wave drag, which is strong due to decreasing air density with altitude and 
significantly inhomogeneous due to randomness of various wave generation and dissipation 
mechanisms. Furthermore, from a global perspective, the adiabatic heating and mechanical 
forcing are balanced in a zonally averaged meridional plane under the quasi-equilibrium 
conditions (Zhu et al. 2001). For example, in the lower stratosphere, because the tropopause is 
much higher (-17 km) in the tropics than in the extratropics (~10 km), an induced thermal 
cooling in the high-latitude lower stratosphere associated with AT t()2 coupled with a mid- 
tropospheric warming in the tropics would enhance a meridional gradient of diabatic heating. 
Such a change in thermal forcing is accompanied by an enhancement in the vertical gradient of 
the wave drag, which is often considered as a dynamical mechanism of driving the strengthening 
of the Brewer-Dobson circulation in the lower stratosphere (Butchart et al. 2006; Garcia and 
Randel 2008; Shepherd and McLandress 2011). 

The equivalent partial change in vertical velocity due to change in COo as shown in Fig. 
6e shows a clear negative correlation to AT t02 shown in Fig. 6a, indicating the fact that a 
decrease in atmospheric temperature can be dynamically associated with an increase in adiabatic 
cooling induced by a strengthening in upward motion. A magnitude of IK in temperature 
decrease due to climate forcing is equivalent to an increase of about 0.025 km day 1 in vertical 
velocity in terms of atmospheric dynamical response. Comparison between Fig. 6d and 6f 
suggests that the linearization from an energy residual to a temperature residual leads to errors of 
less than 10% in partial temperature changes. This can also be considered as a measure of errors 
in converting the generic additive relation (3) for energy differences to the MCFRAM additive 
relation (10) for temperature changes. 

Though AT" 0 " - ™* could be very large locally, its global average in the middle atmosphere 
should be much smaller than its typical local values. This is mainly due to the fact that the 
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globally averaged vertical velocity at a given pressure level should nearly vanish (Olaguer et al. 
1992), and the main role of propagating waves is to redistribute rather than generate the 
momentum and heat (e.g., Zhu et al. 2008, 2010). It is only the eddy diffusion generated by wave 
breaking and molecular viscosity that will be able to produce a globally averaged heating or 
cooling rate difference. In Fig. 7a, we plot the globally averaged partial temperature changes as 
shown in Fig. 6a-d together with the sum of the three components, which gives a very good 
approximation of AT"" 7 in the middle atmosphere. The figure shows that AT'"''" gradually 
increases from 0 near 22 km to IK near 30 km. It remains to be ~1K in the region of 30-70 km. 
The difference between AT'"'" / derived from the direct measurements by SABER and AT'""" is 
At non-md j^ s mean j s shown in Fig. 7b. Figure 7b confirms our conjecture that the 

globally averaged AT""'' is a small difference between globally averaged AT' 0 '" 7 and AT"" 7 in 
most of the middle atmosphere although locally AT" 0 " - ™"* is noticeably greater than either AT'"'" 7 
or AT'""" . Physically, Fig. 7b also means that the globally averaged climate change in the middle 
atmosphere is thermally driven below ~70 km where the vertical eddy transport due to wave 
breaking is expected to be small. An increase in CO 2 concentration coupled with a decrease in 
solar radiation reduces the net radiative heating rate, which cools the global atmosphere. It 
should be pointed out that this is not an obvious scenario among several alternatives. For 
example, an atmosphere can be adiabatically cooled globally at a certain altitude range by a 
systematic upward motion driven by a radiative heating ( e.g ., Zhu et al. 2014). Near and above 
the mesopause region, globally averaged AT""'' "" 7 is no longer small but of the same order of 
magnitude as AT""" 1 or AT"" 7 . This is mainly because the gravity wave breaking in the upper 
mesosphere induces eddy diffusion that irreversibly transports and distributes tracers including 
the potential temperature associated with atmospheric energy. 

It is worth pointing out that the results shown in Fig. 7 independently verify the SABER 
measurements of T and O3 and the accuracy of the JHU/APL radiation algorithm for the middle 
atmosphere. One common way of verifying measurements and testing radiation algorithms is to 
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evaluate the global radiative balance (Kiehl and Solomon 1986; Olaguer et al. 1992). A good 
radiation algorithm requires the globally averaged net radiative heating rate to be much smaller 
than the typical values of the localized net radiative heating rate. A more stringent requirement 
for a good algorithm is to further have a greater sensitivity of heating or cooling rate with respect 
to variations in radiation parameters while still preserving the property of its globally averaged 
net radiative heating rate close to zero. Note that AT"" / and AT''"" md are closely related to the 
difference of the net radiative heating rate and the vertical velocity between two slightly different 
equilibrium states, respectively. The result shown in Fig. 7a suggests that the JHU/APL radiation 
algorithm is sensitive to variations in CO 2 , O3 and F 10.7 and yet the globally averaged AT""” md 
shown in Fig. 7b remains small, as expected for thermally driven global change, based on the 
premise that the SABER measured T and O3 fields are accurate as well. 

4. Application of MCFRAM to GEOSCCM output fields 

Similar to the SABER measurements, we apply the MCFRAM analysis to two- 
dimensional zonal mean fields derived from the GEOSCCM. The 3D GEOSCCM uses the 
GEOS-5 atmospheric general circulation model (Rienecker et al. 2008) in its forecast-model 
component, coupled with the stratospheric chemical solver developed as a part of the GSFC 3D 
chemical-transport model (Douglass et al. 1996; Pawson et al. 2008). With respect to Rienecker 
et al. (2008) this version of GEOSCCM also includes a treatment of stratospheric aerosol (Aquila 
et al. 2012; 2013) and a mechanism to generate the QBO using a gravity wave drag 
parameterization (Molod et al. 2012). The GEOSCCM traditionally uses a fixed input solar 
spectrum, representative of mean solar cycle conditions, and has in fact been used as a no-solar 
cycle reference model in past CCM intercomparisons (Austin et al., 2008). For use in MCFRAM 
the GEOSCCM has been modified to include a solar cycle through the development of new 
atmospheric heating and photolysis code (Swartz et al. 2012). 

In general, the saved fields of GEOSCCM or any other CCMs are not specifically 
designed for directly performing a full MCFRAM analysis. Additional processing to some of the 
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output fields is needed in order to produce a set of appropriate input fields for MCFRAM 
analysis. One potentially important input parameter as shown in Eq. (7) or (9) is the effective 
albedo of surface and the lower atmosphere ( co 0 ) that radiatively couples the middle atmosphere 
with the troposphere and surface. co 0 is not saved in GEOSCCM as an output field in the 
simulations used. We therefore use the saved field “TOA net downward shortwave flux” Frsr 
( W m 2 ) to derive co 0 . Frsr is related to co 0 by the following relationship 

Frsr — S TOA (1 — ) — — ) ’ ( 16 ) 


where S T oa is the TOA downward solar flux with So and /u being the solar constant (= 1366 
W nT 2 , Liou 2002) and cosine of the solar zenith angle, respectively. For a given zonal mean 
F rsr , the diumally averaged S T0A can be calculated by (Cogley and Borucki 1976, Zhu 1994) 


• max(0,zi+5) 


/ud /u 


S T0A = S 0 fi = f ” , 


(17) 


where A = sin ^ sin A, B = cos (j) cos 5 , (j) is the latitude, and 8 is the solar declination angle. We 
finally get the co 0 for applying MCFRAM to a zonal mean field 

( 18 ) 

^TOA 

Note that when A + B < 0 the sun does not rise and S TOA = 0 and co 0 can be any value. Under 
such a circumstance, we set the variation of co 0 between the two states 1 and 2 to be zero. When 
A-B > 0 then the sun does not set and the lower limit of the integration in Eq. (17) is set to 
A-B . Since the upper boundary of the current GEOSCCM is below the mesopause, where the 
effect of O 2 variation is negligible in energy budget, we will neglect AT 02 in this paper. 

Another issue in implementing MCFRAM analysis based on model output fields is that 
most CCMs such as the GEOSCCM only save separately the total solar heating and infrared 
cooling rates but not the individual components contributed by different absorbers and solar flux 
variations. We have already pointed out previously through Eq. (6) that the separation of cooling 
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rate components in the middle atmosphere is very simple because the effect of line overlapping 
is negligible. On the other hand, Eq. (7) suggests that it requires a significant overhaul to the 
online radiation code in any CCMs in order to derive and save heating rate contributions by 
different components mainly because of the nonlinear effect between solar flux and absorber. 
One way to get around the whole issue is to calculate all the radiative heating and cooling 
perturbation terms offline and introduce two error terms to the basic MCFRAM Eqs. (10) and 
(1 1) (Taylor et al. 2013; Sejas et al. 2014): 

AT lotoi _ ^ at ( ' !) - AT e,rl - AT" 7 ” 2 . (19) 

n 

Here, the two partial temperature changes due to radiation errors are calculated based on 
GEOSCCM-saved total solar heating and infrared cooling rates together with the offline 
radiation algorithm: 

AT OTl =Z-AS OT and AT c,r2 = Z-(-AR" T ) , (20a, b) 

where AS' :7T and AR' rr are respectively the changes in total radiative heating and cooling rates 
between two states 1 and 2 derived from the offline and GEOSCCM online radiation algorithms 

AS OT = AS <# -AS ecm and AR OT = AR o// - AR ccm . (21a,b) 

It has been suggested that the error terms are mostly contributed from the different averaging 
procedures between the online and offline calculations (Taylor et al. 2013; Sejas et al. 2014). 
Additional errors will also contribute to partial temperature changes due to radiation errors when 
different radiation algorithms are adopted for the online and offline radiative heating and cooling 
rate calculations. The error introduced by inferring AQ from radiative forcing evaluated from 
the offline radiative transfer model calculations can be estimated and analyzed by a comparison 
with that derived directly from CCM outputs saved during runtime (Sejas et al. 2014). 

In this paper, we choose the same output time periods of 2002-2003 (near solar 
maximum) and 2008-2009 (solar minimum) from one GEOSCCM simulation as those for 
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SABER measurements used in the last section to perform the MCFRAM analysis. In Fig. 8, we 
show the variation in effective albedo of the surface and lower atmosphere scaled by the 
diurnally averaged solar radiation (S 0 ju)Aco 0 as a function of month and latitude over the 24- 


averaged solar flux ( S 0 / 4 ) and is one order of magnitude greater than the variation in solar 
constant for the 11 -year solar cycle (Fean 1991). The figure shows significant geographic and 
transient variations with peak values appearing near equatorial and summer polar areas where the 
maximum mean solar fluxes are deposited. Climate change or the system’s feedback response is 
often associated with a radiative forcing scaled by the changes in the total radiation flux. Since 
the energy deposition in the atmosphere at different wavelengths varies drastically with spatial 
and temporal distributions of absorbers, the change in the input solar energy may not be able to 
fully represent how the system responds. On the other hand, the MCFRAM analysis based on 
Eqs. (10)-(1 1) together with Table 1 provides us with a complete view of the system response in 
the same variable and units under an observational constraint of the measured total temperature 
change (AT ro ' fl '). 

In Fig. 9, we show all the partial temperature change components in Table 1 in the middle 
atmosphere below 70 km that can be directly calculated based on GEOSCCM output fields and 
the offline JHU/APF radiation algorithm. Panels (a)-(j) correspond to the first 10 rows in Table 1 
plus Eqs. (20a, b) but excluding AT 02 , which is negligible below the mesopause. Panels (k) and 
(1) are respectively the partial temperature changes due to all radiative processes ( AT"" 7 ) with the 
offline radiation algorithm only ( AT'^ n( , ) and that including the online correction terms 


month period. The figure shows a typical variation of ~5 W m 2 that is about 1% of the globally 




( 22 ) 


Panel (m) sums all the AT (,!) components that can be directly calculated 


AT' 


''Sum 

online 


= at;1 +at w *+at®. 


online 


(23) 
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Panels (n) and (o) are the residual partial temperature changes corresponding to the online 
versions of AT etW> ' and AT""" rad defined in Table 1 , respectively. 

We first note that the overall patterns and magnitudes of the partial temperature changes 
AT C02 and AT F107 (Figs. 9a and 9c) that are primarily induced by the variations of the external 
forcing are nearly identical to those derived by SABER measurements in the common domain 
(Figs. 6 a and 6 c). However, GEOSCCM shows an additional strengthening in partial temperature 
change associated with the CO? cooling in the available high latitude and polar regions, 
especially in the southern hemisphere mesosphere where the coldest temperature often occurs 
near the summer mesopause (Lubken 1999; Lubken et al. 1999). This is caused by a non- 
localized heat exchange between the warmer stratopause and colder mesopause when the CO 2 15 
pm band transmission behaves transparently and the summer mesopause receives net radiative 
heating from the stratopause (Zhu et al. 1992). An increase in CO 2 concentration increases the 
atmospheric opacity that leads to a reduction in summer mesopause net heating rate. 
Furthermore, there exists a local maximum in CO 2 15 pm cooling rate near the winter polar 
mesopause due to the combination of local thermodynamic equilibrium conditions, a near 
uniform temperature, and a near transparent emission to space (Zhu 1994). This too contributes 
to the strengthening in AT to: in the high-latitude and polar mesosphere regions. There exists a 
significant difference in AT 03 between GEOSCCM fields (Fig. 9b) and SABER measurements 
(Fig. 6 b). This is not surprising because the middle atmosphere O3 and its variability are very 
sensitive to a strong nonlinear coupling between photochemistry and dynamics. For example, the 
largely off-set peaks in AT 03 in the equatorial lower stratosphere may well reflect the degree of 
fidelity of GEOSCCM simulation to the equatorial quasi-biennial oscillation phenomenon. 

The partial temperature change AT H2 ° as shown in Fig. 9d makes a much smaller 
contribution and is negative in the low latitude but positive in part of the midlatitude in the 
middle atmosphere. Middle atmosphere H 2 O may increase with time as a result of increasing 
CH 4 in the troposphere (Zhu et al. 1999). Its decadal change could also be well correlated to the 
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equatorial sea surface temperature that largely determines the coldness of tropopause to limit the 
direct entry of H 2 O into the stratosphere (Solomon et al. 2010). The existence of large regions of 
both positive and negative AT H2 ° in the middle atmosphere is an indication that both processes 
play important roles in determining HiO concentration in the time period of 2002-2009. The 
contribution by the lower atmosphere effective albedo AT'" 0 is even smaller than AT H2 ° by 
nearly one order of magnitude and changes are largely confined in the stratosphere (Fig. 9e). 
Comparison between Fig. 8 and Figs. 9a-e gives us one example that the MCFRAM with its key 
additive property provides us with a more direct and quantitative insight into the relative 
importance of different factors of climate forcing and feedback processes when they are 
constrained under the same scale with the same units. Panels (f) and (g) represent the part of the 
dynamical effects on the atmospheric thermal response that can be easily evaluated based on the 
available model output fields. As we have already conjectured in discussing the MCFRAM 
applications to the SABER measurements, the directly calculable components of AT"'"' r ‘"' are 
overwhelmingly large with the majority of the contributions coming from AT 1 '* . We note that 
the peak values of AT" occur near polar areas and arctic and antarctic circles, where the 
spherical geometry may lead to unusually large variability in solar heating rate or flow 
divergence, both in the real atmosphere and in numerical models. 

Panels (h) and (i) in Fig. 9 show the partial temperature changes, AT"” 1 and AT" ,r2 , due 
to differences in heating and cooling rates between the offline and online calculations, 
respectively. The figures show that the differences are small in most regions of the middle 
atmosphere except AT 6 "’ 1 near the low latitude upper boundary and AT 6 " 2 near the polar area. 
We note that the heating rate difference associated with AT 6 " 1 by the solar radiation near 
model’s upper boundary is sensitive to the shielding effect of the solar flux by the absorber 
column above the upper boundary. Furthermore, the sensitivity decreases with increasing latitude 
as the slant path also increases. The high latitude cooling rate difference associated with AT 6 " 2 
is likely sensitive to the non-localized heat exchange when the vertical temperature gradient is 
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large. Specifically, the heat exchange by the CCb 1 5 -pm band becomes transparent above the 
stratopause for a uniform CO2 mixing ratio distribution whereas the O3 9.6 pm band emission 
could be transparent in the entire atmosphere in regions where the O3 concentration is low (Zhu 
et al. 1991, 1992, Zhu 1994). The issue of AT e ' rl and AT" 7 ’ 2 will be further pursued in the next 
stage of investigation. 

Figure 9j shows the measurable temperature difference ( AT' , "" / ) between the two 
equilibrium states 1 and 2, which is the difference of the model output temperature fields. 
Comparing with Fig. 5c, we note that both the modeled and SABER measured AT' , "" / show 
positive-negative paired peaks of the same magnitudes near 50° latitudes and equatorial lower 
stratosphere. We note that AT'"'" 7 provides an observational constraint and a unified or a 
standard scale to all the other sensitivity responses in the MCFRAM analysis. Recall that 
MCFRAM (or CFRAM) was developed from the energy budget equation (1) or (2). In addition 
to changes in the energy budget due to all the parameter variations shown in panels (a)-(g), the 
most prominent and well behaved one is the change in cooling rate and the diffusive heat 
exchange caused by the variation of atmospheric temperature T. The well-behaved nature of 
A(R-Q mo/ ) with respect to T in Eq. (4) makes the generalized damping matrix A always 
invertible. Furthermore, the temperature T is a directly measurable and most common variable. 
These two features can be considered the underpinning for MCFRAM that exclusively separates 
the temperature component of variation in the energy budget from all the other components in 
Eq. (4) and set it to be a standard scale to be compared to all the other feedback responses. 
Panels (k) and (1) in Fig. 9 show the partial temperature changes due to radiative processes based 
on offline and online radiation algorithms, AT^ and AT™^ ( = AT^ - AT OTl - AT^ 2 ), 
respectively. We note that the magnitude of AT"" 7 increases with altitude, which is consistent 
with that of the measurable AT'""' / . However, there exist significant differences in spatial 
structure contributed from the partial temperature changes due to non-radiative processes 
Ay non-rad gj ncc ^non-rad j g dependent on atmospheric motion that is strongly nonlinear and 
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contains many different scales, we expect the magnitude AT'"' 1 rad to be reduced when an 
ensemble average is taken for the MCFRAM analysis to the output fields from many different 
runs of GEOSCCM in our future investigations. 

Among all the components AT (,,) listed in Table 1 the biggest component that can be 
directly calculated as shown in Fig. 9f is AT" . Its typical localized value is nearly an order of 
magnitude greater than temperature changes AT 1011 ' 1 or AT' m/ as shown in Fig. 9j-91. As a result, 
if we sum up all the terms in Table 1 that can be directly calculated, AT 11 M as shown in Fig. 9m 
in its online version, then its overall spatial distribution will be dominated by that of AT" - * . The 
last two panels (n) and (o) in Fig. 9 correspond to two last partial temperature changes in Table 
1, AT ai * and AT'“' ! rad , calculated by the residual method. AT 0 **’’ is due to dynamical heating 
contributed by the unresolved eddies and horizontal winds. AT” _,W is /\T Kddv plus the partial 
temperature changes associated with the adiabatic cooling due to vertical motion that can be 
calculated based on the column profiles. We see again that AT' yW> largely cancels AT"* 
contained in AT 11 S) mostly due to the energy perturbation associated with the horizontal 
motions. Since the localized partial temperature changes due to radiative processes AT'"' / shown 
in panel (1) is smaller than the local values of AT' oto/ shown in panel (j), the overall magnitude 
and spatial structure of AT '“" ~ md as shown in the last panel (o) is similar to those of AT'"' a/ , 
indicating dynamical processes dominate the local structure of the total partial temperature 
change. This is also consistent with our previous analysis to SABER measurements where Fig. 
5c and Fig. 6d show large similarities in their overall magnitude and spatial structure. 

In Fig. 10, we show the globally averaged partial temperature changes presented in Fig. 
9. Several major features are consistent with those derived from the SABER measurements as 
shown in Fig. 7: (i) AT F107 makes the largest contribution above —40 km, (ii) AT C ° 2 is negative 
at all altitudes whereas AT 03 is positive in some part of the altitude range, (iii) the globally 
averaged AT" due to atmospheric circulation makes a small contribution to the global mean 
climate change in the middle atmosphere below 70 km, which is primarily driven by radiative 
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processes. In addition to AT" , the three more terms AT H2 °, AT® and AT® 0 have been directly 
calculated based on the GEOSCCM output fields. Since the globally averaged values of these 
terms are all smaller than KT rad , our major results derived from the SABER measured T and O3 
fields remain valid. We have already mentioned that it is generally unavoidable to adopt an 
offline radiation algorithm to perfonn the MCFRAM analysis. We note that AT 6 ” 2 due to 
cooling rate difference is negligibly small, and AT'' ,rl increases rapidly near the model’s upper 
boundary, which in turn leads to a large deviation of AT"" 7 from AT'"'" 7 near the upper 
boundary. We note that all the model fields in GEOSCCM, including the atmospheric 
temperature, have been integrated subject to the influence of a set of prescribed boundary 
conditions. On the other hand, the magnitude of AT"" 7 derived from the SABER measurements 
as shown in Fig. 7 does not systematically increase with the altitude below 80 km, indicating the 
effect of boundary condition on the heating rate calculations for GEOSCCM fields. We will 
pursue this issue in our follow-up investigations. 

5. Summary 

In this study, we have extended the Climate Feedback-Response Analysis Method 
(CFRAM) for the coupled atmosphere-surface system to the middle atmosphere. The Middle 
atmosphere CFRAM (MCFRAM) is built upon the atmospheric energy equation per unit mass 
with radiative heating and cooling rates as its major thermal energy sources. In addition, 
molecular thermal conduction is added to the energy equation when the upper boundary is 
extended beyond the mesopause. MCFRAM preserves the unique feature of an additive property 
for the original CFRAM in which the sum of all partial temperature changes equals the observed 
temperature change. By introducing the generalized damping (A) and relaxation (Z) matrices to 
the basic MCFRAM equation, the relationship between the fundamental quantity of the partial 
temperature change (AT (n) ) and its physical cause of energy perturbation (AF (,,) ) is 
quantitatively clarified by the well-documented theory of radiative damping of thennal 
disturbance in the middle atmosphere. Specifically, we show that A serves as a filter that 
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smoothes the small-scale structure in AF ( ”\ In addition, it is shown that for a given energy 
perturbation the maximum response in temperature change occurs when the energy perturbation 
is located at the place where the cooling rate of the mean state reaches it minimum value. 

The newly developed MCFRAM is applied to two sets of two-dimensional data. One is 
the zonal mean T and O 3 fields in the middle atmosphere derived from SABER measurements in 
the low and midlatitudes averaged over yaw cycles. The other is the zonal mean fields saved 
from GEOSCCM simulations. It is found that the spatial structure of the temperature responses 
to variations of CO 2 , O 3 and solar flux are different. AT C02 closely follows temperature 
distribution in most of the middle atmosphere because the cool-to-space approximation is valid 
for an atmosphere with uniformly distributed CO 2 mixing ratio. Both the solar radiation heating 
and 9.6-pm band cooling by O 3 affect AT 03 in about the same order of magnitude, both 
processes strongly influenced by O 3 distribution. AT FI07 monotonically increases with altitude 
due to the fact that the solar UV fluxes of greater variations at shorter wavelengths are generally 
absorbed at higher altitudes. The two periods used to derive the statistical equilibrium states are 
2002-2003 and 2008-2009, corresponding to near solar maximum and solar minimum, 
respectively. The CO 2 mixing ratio between these two periods increases from -374.7 ppmv to 
~386.3 ppmv. It is consistently found by both datasets that for a half cycle span of the 11-year 
solar cycle the largest component of the partial temperature changes (AT (,,) ) in the middle 
atmosphere is the one due to the variation of the input solar flux (AT F107 ). The effect of 
increasing CO 2 always cools the middle atmosphere with time (AT to 2 <0). On the other hand, 
depending on the relative importance of O 3 heating and cooling rates, AT 03 could be either 
positive or negative. The MCFRAM analysis to GEOSCCM fields suggests that AT H2 ° makes a 
minor contribution to the total temperature change observed from the atmosphere ( AT""" / ). The 
partial temperature change due to the variation of the effective albedo of the surface and lower 
atmosphere to the solar radiation (AT“°) is negligibly small in comparison with those by other 
factors. 
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Because of the lack of all the required parameters in the input datasets, the partial 
temperature change due to non-radiative processes ( ^non-rad ^ 0 q cn necc i s to be evaluated by 
either an energy or a temperature residual approach. Such an approach is well-founded due to the 
existence of the additive property for the generic energy equation (3) or the basic MCFRAM 
equation (10) to temperature changes. AT"°"~™ d for the SABER measurements includes all 
dynamical effects whereas three individual components in AT"°' , ~''“ / can be evaluated separately 
based on the GEOSCCM model outputs. In both cases, the typical magnitude of AT Mj is 
significantly greater than any component consisting of partial temperature changes due to 
radiation processes ( AT'" £/ ). However, the global average of AT""" ™ ,/ is much smaller than that 
of AT"" / below -70 km, indicating the lack of vertical transport of energy by eddies or by global 
mean vertical velocity. Physically, this means that the globally averaged climate change in the 
middle atmosphere below -70 km is thermally driven. This also means that the globally 
averaged partial temperature change due to all radiative processes is approximately equal to the 
observed temperature change. It ranges from -0.5 K near 25 km to -1.0 K near 70 km from the 
near solar maximum to the solar minimum. 
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Z-(AS 03 -AR 03 ) 
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z-as F107 

AT (3> due to change in downward solar radiation at TOA 
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Z-(-AR H2 °) 

AT' 41 due to changes in H^O 
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Z-AS 02 

AT 1 '' 1 due to changes in 0 2 

AT® 0 

Z • AS® 0 

AT (6) due to changes in troposphere albedo to the solar radiation 
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Z • (0Aw*) 

AT (7) due to changes in the resolved vertical velocity 
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Z • (A0)w* 

AT* 81 due to changes in the static stability 
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Z • (A"Q erf<<v ) 

AT <g ’ due to changes in the un-resolved eddies 

non— rad 
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FIGURE CAPTIONS 


Figure 1. Global mean temperature and ozone profiles in the middle atmosphere derived from 
TIMED/SABER measurements in the low and mid-latitudes over a 12-year period of 2002-2013. 
The TIMED/SABER measurements in the middle atmosphere are merged to the US Standard 
Atmosphere in the troposphere. The vertical resolution is about 0.7 km. 

Figure 2. Selected vertical eigenmodes of the generalized damping matrix A calculated from T 
and O 3 shown in Fig. 1 based on the JHU/APL radiation algorithm. The CO 2 volume mixing 
ratio is set at 2005 level of 380 ppmv. The unit of the eigenvalues shown in the figure boxes is 
day -1 . 

Figure 3. A quantitative relationship between the generalized damping rate and the vertical 
wavenumber at which the poser spectral density is maximally peaked. Also shown in the figure 
are analytic fits of radiative damping given by (Zhu 1993) and a fit for diffusive damping 
proportional to the square of the vertical wavenumber. 

Figure 4. Linear temperature responses to three energy perturbations caused by changing three 
atmospheric parameters (i) CO2 volume mixing ratio is doubled from 380 ppmv to 760 ppmv, 
(ii) O3 mixing ratio is uniformly reduced by 50%, and (iii) solar flux index F10.7 is increased 
from 60 to 260 (in units of 10 ~ 22 W m 2 Hz ! ). 

Figure 5. Zonal mean T and O 3 fields averaged over a 12-year period of 2002-2013 and their 
differences between two equilibrium states covering the periods of 2002-2003 and 2008-2009, 
respectively. 

Figure 6 . Two dimensional distributions of partial temperature changes between two time 
periods of 2002-2003 and 2008-2009 due to variations in (a) CO 2 , (b) O 3 , (c) F 10 . 7 , and (d) 
atmospheric circulation, respectively, (e) Equivalent partial change in vertical velocity due to 
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change in CO?, (f) Error in partial temperature change of non-radiative processes due to 
linearization approximation. 

Figure 7. Globally averaged partial temperature changes AT C02 (dashed line with squares), 
AT 03 (dashed line with triangles), AT F107 (dashed line with circles) and their sum approximately 
representing KT vad (solid line with solid circles), (b) Globally averaged total temperature change 
AT"" fl/ (dashed line) and partial temperature changes AT’™' (solid line with solid circles) and 
Ay non-rad ( so |j c | | me w j th diamonds). 

Figure 8. Changes in effective albedo of the surface and lower atmosphere scaled by the 
diurnally averaged solar radiation between 2002-2003 and 2008-2009. 

Figure 9. Partial temperature changes caused by various energy perturbation components in the 
middle atmosphere. The unit in scale bars of all panels is K. 

Figure 10. Globally averaged partial temperature changes shown in Fig. 9. 
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(a) Global mean T profile measured by SABER 



C > 3 volume mixing ratio 


(b) Global mean O 3 profile measured by SABER 



962 Figure 1. Global mean temperature and ozone profiles in the middle atmosphere 

963 derived from TIMED/SABER measurements in the low and mid-latitudes over a 12- 

964 year period of 2002-2013. The TIMED/SABER measurements in the middle 

965 atmosphere are merged to the US Standard Atmosphere in the troposphere. The vertical 

966 resolution is about 0.7 km. 
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974 Figure 2. Selected vertical eigenmodes of the generalized damping matrix A calculated 

975 from T and O 3 shown in Fig. 1 based on the JHU/APL radiation algorithm. The CO 2 

976 volume mixing ratio is set at 2005 level of 380 ppmv. The unit of the eigenvalues 

977 shown in the figure boxes is day -1 . 
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Damping rate versus vertical wavenumber 



Vertical wavenumber (radian km" 1 ) 


981 

982 Figure 3. A quantitative relationship between the generalized damping rate and the 

983 vertical wavenumber at which the poser spectral density is maximally peaked. Also 

984 shown in the figure are analytic fits of radiative damping given by (Zhu 1993) and a fit 

985 for diffusive damping proportional to the square of the vertical wavenumber. 
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Figure 4. Linear temperature responses to three energy perturbations caused by 
changing three atmospheric parameters (i) CO 2 volume mixing ratio is doubled from 
380 ppmv to 760 ppmv, (ii) O 3 mixing ratio is uniformly reduced by 50%, and (iii) 
solar flux index F 10.7 is increased from 60 to 260 (in units of 10 22 W nT 2 Flz 1 ). 
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Figure 5. Zonal mean T and O 3 fields averaged over a 12-year period of 2002-2013 and 
their differences between two equilibrium states covering the periods of 2002-2003 and 
2008-2009, respectively. 
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(b) Mean AT due to AfOJ 
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Figure 6 . Two dimensional distributions of partial temperature changes between two 
time periods of 2002-2003 and 2008-2009 due to variations in (a) CO 2 , (b) O 3 , (c) 
F10.7, and (d) atmospheric circulation, respectively, (e) Equivalent partial change in 
vertical velocity due to change in CO 2 . (f) Error in partial temperature change of non- 
radiative processes due to linearization approximation. 
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Figure 7. (a) Globally averaged partial temperature changes AT C ° 2 (dashed line with 
squares), AT 03 (dashed line with triangles), AT FI07 (dashed line with circles) and their 
sum approximately representing AT™' 7 (solid line with solid circles), (b) Globally 
averaged total temperature change AT"' to/ (dashed line) and partial temperature changes 
AT"" 7 (solid line with solid circles) and AT'"”' rad (solid line with diamonds). 
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Figure 8. Changes in effective albedo of the surface and lower atmosphere scaled by the 
diurnally averaged solar radiation between 2002-2003 and 2008-2009. 
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Figure 9. Partial temperature changes caused by various energy perturbation 
components in the middle atmosphere. The unit in scale bars of all panels is K. 
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Figure 10. Globally averaged partial temperature changes shown in Fig. 9. 
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